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The classical Heisenberg antiferromagnet on a triangular lattice with the single-ion anisotropy of 
the easy-axis type is studied theoretically. The phase diagram in an external magnetic field is con- 
structed from the mean-field analysis. Three successive Berezinskii-Kosterlitz-Thouless transitions 
are found by Monte Carlo simulations in zero field. Two upper transitions are related to the break- 
ing of the discrete Ze-symmetry, while the lowest transition is associated with a quasi-long-range 
ordering of transverse components. The intermediate coUinear phase between the first and second 
transition is the critical phase predicted by J. V. Jose et al. [Phys. Rev. B 16, 1217 (1977)]. 
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I. INTRODUCTION 

Frustrated magnetic systems have been a stimulating 
research topic over several decades. Their diverse prop- 
erties, highly degenerate ground states, non-coUinear or- 
dering, novel phase transitions,i offer a playground to 
investigate fundamental physical questions going far be- 
yond magnetism itself. One of the specific subjects in this 
field is the interplay of geometrical frustration and mag- 
netic anisotropies. The prominent example is provided by 
the rare-earth pyrochlore materials with Ising-type mag- 
netic moments. Contrary to naive expectations, these 
magnetic systems remain non-frustrated for an antiferro- 
magnetic nearest-neighbor coupling, but develop highly 
frustrated spin-ice states for the case of a ferromagnetic 
exchange between spins i^i^ii^ 

In the present work we investigate the nearest-neighbor 
Heisenberg antiferromagnet on a triangular lattice with 
the single-ion anisotropy of the easy-axis type: 

n = jY,s.-s,-Dj2istf (1) 

{ij) i 

Such a Hamiltonian is believed to describe quasi two- 
dimensional (2D) magnetic materials VCI2 (Ref. [Sj) and 
LiCr02 (Ref. i). A similar model with the XX Z 
anisotropy has been previously studied by a number of 
authors. ^'^^iS In real magnetic materials with 5* > i the 
single-ion anisotropy being the first-order relativistic ef- 
fect is usually more significant than the anisotropic ex- 
change, which is generally of the second-order in the spin- 
orbital coupling.— Besides, as we shall see later, the two 
types of anisotropy lead to different sequences of finite- 
temperature phase transitions. 

Ordered states of the anisotropic triangular antiferro- 
magnet ([TJ are characterized by a nonzero static magne- 
tization: 

(S,) = li cos(Q • r,) + h sin(Q • r,) + m . (2) 

with the ordering wave vector Q — (47r/3, 0). At zero 
temperature the Heisenberg triangular-lattice antiferro- 
magnet orders in a three-sublattice 120° spin structure. 




FIG. 1: (Color online) Possible three-sublattice planar con- 
figurations of the easy-axis triangular antiferromagnet. The 
direction of the easy-axis is shown by z. The nonzero compo- 
nents of the order parameter ((2]) are indicated for each con- 
figuration. 



Such a noncoUinear magnetic ordering is described by 
a pair of orthogonal antiferromagnetic vectors: h _L I2, 
|li| — II2I, and m = 0. In accordance with the Mermin- 
Wagner theorem there is no symmetry breaking tran- 
sition at any finite temperature. Still a weak topo- 
logical transition related to proliferation of Z2-vortices 
may occur for this model at T/J - O.SJii^ii^iiiii^ For 
the easy-plane anisotropy, Z? < in Eq. the spin 
plane of the ordered 120° structure is fixed to the x-y 
plane. In this case two finite temperature transitions take 
place: the Ising-type transition related to the chiral sym- 
metry breaking and the Bcrczinskii-Kostcrlitz-Thouless 
(BKT) transition associated with the vortex-antivortex 
unbinding. 

The easy-axis anisotropy, > 0, orients the spin plane 
perpendicular to the x-y crystallographic plane and si- 
multaneously distorts the spin structure. Finding direc- 
tions and magnitudes of li and I2 becomes a nontrivial 
problem in this case. Possible spin structures correspond- 
ing to the ordering wave vector Q are presented in Fig. [T] 
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They have been obtained by a symmetry analysis and 
are confirmed by the mean-field calculations described in 
the next section. Some of these states, Figs. HJa), (c), 
and (d), have a finite uniform magnetization m along z, 
which is, however, a secondary order parameter and not 
indicated for that reason in the Figure. 

In order to elucidate symmetries of different phases, we 
note that a simple translation Ta (r.; — > r^+a) transforms 
the antiferromagnetic order parameter according to 

Ta[(ll+^l2)] = ih+th) e-'''^-" , (3) 

where the phase factor can take only three different val- 
ues: Q • a = 0, ±27r/3. Hence, besides the group of 
continuous rotations about the z-axis the magnetic struc- 
ture has an inherent discrete symmetry Z3. Such an ad- 
ditional symmetry corresponds to permutations of three 
sublattices. In zero magnetic field the time-reversal sym- 
metry implies invariance with respect to li —li, which 
enlarges Z3 to Zg. The total symmetry group is, there- 
fore, 

G = 5i «) Ze , (4) 

see also a similar discussion in Ref. j^. The coUinear 
phases shown in Figs. [Ija)-(c) preserve the axial sym- 
metry 5*1 but break in different ways the discrete sym- 
metry group Zg. In terms of the order parameter angle 
if defined by 

hz^lcosip, hz^lsinip, (5) 

the state in Fig. [TJa) corresponds to commensurate val- 
ues If — 2fc7r/6 with an integer k, whereas the configu- 
ration in Fig. [TJb) has (p = {2k + l)7r/6. The third type 
of a coUinear state is described by an arbitrary angle ip 
and is shown schematically in Fig. [He). In such a state 
the phase ip remains unlocked and the sine and cosine 
harmonic ([S]) coexist with an arbitrary ratio. 

For large enough values of D/J > dc = 1.5 the mag- 
netic anisotropy induces a highly degenerate coUinear 
Ising state at zero temperature. Quantum fiuctuations 
can lead, then, to interesting zero- and finite-temperature 
phases Ji'-lS Here, we investigate an antiferromagnet with 
a moderate-strength anisotropy Q < D / J < dc, which is 
frequently found among experimental systems, and con- 
sider the finite-temperature properties of the model ([1]). 
For simplicity, we neglect quantum effects and study the 
classical spin model. 

Layered easy-axis triangular antiferromagnets with a 
significant interplane coupling exhibit two second-order 
transitions with an intermediate coUinear /i-phase shown 
in Fig. [ija).''^^ In contrast, we show in the present work 
that a purely 2D system ([T|) shows three consecutive 
BKT-type transitions. In the first part, Sec. II, we inves- 
tigate the mean-field phase diagram of the model H]) at 
zero and at finite magnetic fields. The mean-field behav- 
ior is expected to be realized in layered triangular anti- 
ferromagnets with weak interplane coupling. The Monte 
Carlo (MC) simulations and the analysis of the zero-field 



behavior of the model ((T)) are presented in the second 
part of our study, Sec. III. 

II. MEAN-FIELD THEORY 

Let us begin with the mean-field analysis of possi- 
ble finite-temperature phases of the model (H)). Specif- 
ically, we use the real-space approach;22i2ii^i2^ gener- 
alizing the previously established technique to systems 
with the single-ion anisotropy. The two standard steps 
of the mean-field approximation include (i) decoupling 
the spin-spin interaction according to 

S, . S, « S, . (S,> + (S,) • S, - (S,) . (S,) , (6) 

with (Si) being the thermal average of an i'^ magnetic 
moment and (ii) rewriting 7i as a sum of single-site 
Hamiltonians 

Hmf - 51 [ ~ i?(^f )' - h, • S,] - J^(S,) • (S,) 

with h, = H - ./^ (Sj) , (7) 

11. n. 

where we have also added a Zeeman magnetic field to 
Eq. JT]). Due to the presence of the single-ion term in 
TYmf, the local magnetization (Si) has to be decomposed 
into components, which are transverse and parallel to the 
anisotropy axis: 

{S,) = {St)i+{St) ^^'"^'^^ . (8) 

K 

Performing integration with respect to x = Sf = cos 9i 
in the expression for the partition function we obtain the 
following mean-field equations for static magnetic mo- 
ments: 

{SD - ^J'dxxe^^'/^e'^>/^ I,{y.), (9) 

where yi = h^y/l — /T and /„(z) is the modified 
Bessel function of the n-th order: 

/„(z) = i / dae^™^" cos" a . 
I" J a 

The system of integral equations ([9]) together with the 
self-consistency condition given by Eq. ([7]) is solved iter- 
atively on finite lattices of N = L x L spins, with peri- 
odic boundary conditions. Once convergence is achieved, 
various physical quantities are calculated including the 
free-energy 

.FMF = -J^(S.)-(Sj)-r^lnZ, , (10) 
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FIG. 2: (Color online) The low-field part of the mean-field 
phase diagram with H || z for a representative value D/ J = 1. 
Spin configurations for each phase are schematically indicated 
by arrows as in Fig. [1] Solid and dashed lines correspond to 
first- and second-order transitions, respectively. 

the internal energy Emf — {'Hmf), and the antiferro- 
magnetic order parameters. By exphcit calculations for 
clusters with 3 < L < 12 at all temperatures and weak 
magnetic fields we have verified stability of the three- 
sublattice structure with Q = (47r/3,0). After that a 
more detailed investigation of the H-T phase diagram 
has been performed with the three-sublattice ansatz. 
Precise location of phase boundaries in Fig. [2] has been 
determined from temperature and field scans for the an- 
tiferromagnetic order parameters indicated in Fig. [1] as 
well as for the uniform magnetization. The behavior of 
the specific heat has been also used to independently ver- 
ify these results. 

At the upper transition Td ~ 1.2 J in zero magnetic 
field only z-components of magnetic moments become 
ordered. In accordance with the Zg symmetry selection 
between various collinear structures is determined by the 
following invariant in the Landau free-energy: 

Aeiilt+tm'^+c.c.] . (11) 

For negative Aq < the pure Zi-state, Fig. [TJa), is en- 
ergetically favored, while Aq > corresponds to the I2- 
state, Fig.[Ub). We have verified the positive sign of Aq 
in our case by a direct analytical expansion of Eqs. ([9]). 
Our numerical results also confirm that the Z2-state is 
stable below Td. Such a partially ordered phase has a 
vanishing moment on one of the antiferromagnetic sub- 
lattices. A similar phase has been discussed in relation 
to the intriguing phase diagram of Gd2Ti207.— Here, we 
provide an example, where a partially ordered phase is 
realized at the mean-field level in a simple spin model. 

The second transition at Tc2 — 0.6J is related to the 
breaking of the rotational symmetry about z-axis. Below 
Tc2 the third previously disordered magnetic sublattice 
becomes ordered with moments oriented within the x-y 
plane. Simultaneously, moments of the other two sublat- 
tices start deviating from z-axis leading to a distorted tri- 
angular structure shown in Fig. [IJe) . This distorted spin 



structure is characterized by I2 || z and li _L I2. When 
temperature is further decreased the coefficient Aq in the 
effective anisotropy term changes sign at Tcs ~ 0.3 and 
one finds a first-order transition into another distorted 
triangular structure shown in Fig. [ijd) with li || z. 

Note, that the related model with the exchange 
anisotropy''''* has Aq ^ in the mean-field approxima- 
tion, which leads to an additional continuous degener- 
acy. As a result, only two finite-temperature transitions 
are found in this case: from the paramagnetic state to 
a degenerate collinear configuration shown in Fig. [Ijc) 
and then to a degenerate distorted 120° configurationiiS 
Sheng and Henley^ have discussed how different types 
of fluctuations, thermal, quantum, or random dilution, 
can induce a finite Aq. For the model with the single- 
ion anisotropy one finds a different interesting possibility: 
the sign of the anisotropic term changes upon lowering 
temperature. 

The two phases in Figs.UJa) and (d) have a nonvanish- 
ing total magnetization . The coupling between ferro- 
and antiferromagnetic components is determined by the 
term 

m''{ll+iqf + C.C. , (12) 

which is invariant under Z3 transformations ([3|). In zero 
magnetic field this yields ~ (Tc — T)^/^ for states with 
li II z. In contrast, states in Fig.[IJb) and[lje) with I2 || z 
have vanishing . This difference is important to under- 
stand the finite-field behavior, see Fig. [2l Magnetic field 
applied parallel to the z-axis favors spin structures with 
a finite magnetization and stabilizes states with If ^ 0, 
which is why the two intermediate low-field phases are 
no longer pure Z| states. This feature is emphasized by 
hatches in Fig. O The coUinear-noncoUinear transitions 
are of the second order, whereas all other transition lines 
are of the first order. In the case of the transition from 
the paramagnetic state in external magnetic field the 
first-order nature of the transition follows from the pres- 
ence of the cubic invariant ([H]), while in other cases the 
above conclusion is a consequence of the group-subgroup 
relation. The transition lines intersect at a multicritical 
point {T*,H*) = (0.6 J, 0.25 J). 

The mean-field phases and the structure of the phase 
diagram at fields larger than H* are similar to the Heisen- 
berg triangular antiferromagnet^- so we do not go into 
further details. We have also checked other moderate val- 
ues of D/J < 1.5 and found precisely the same structure 
of stable phases with triple transitions in zero magnetic 
field. As we shall see in the next section, the true ther- 
modynamic phases determined by Monte Carlo simula- 
tions of the model ([T]) differ from the mean-field solutions, 
which is often the case in 2D. Still, the mean-field pic- 
ture is expected to be qualitatively correct for 3D layered 
triangular antiferromagnets. By including a ferro- or an- 
tiferromagnetic interlayer coupling J' in the mean-field 
equations ([7]) and ^ we have verified that the predicted 
sequence of finite-temperature transitions remains valid 
up to I J'/>^l ^ 0-6. For larger values of | J'/>^l we find a 
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FIG. 3: Schematic zero-field phase diagram of the two- 
dimensional triangular antiferromagnet with easy-axis single- 
ion anisotropy. The arrow labeling of phases is the same as 
in Fig. 1. 

double transition with an intermediate h collinear phase 
similar to the previously studied case of very strong J'J^ 

III. MONTE CARLO SIMULATION 

In uniaxial magnetic systems, transverse and longitu- 
dinal spin components order at different temperatures 
as they belong to different irreducible representations. 
For the triangular antiferromagnet with the easy-axis 
anisotropy, the highest transition should be related to the 
sole breaking of Zg symmetry. Such a discrete symmetry 
breaking may lead to a phase with a true long-range or- 
dering at low temperatures even in 2D. The case of a 2D 
system with the general Zp symmetry has been consid- 
ered in the seminal work of Jose and co-workers.^^ The 
precise nature and sequence of finite-temperature transi- 
tions depend on the number p of "clock states." Jose et 
al. have predicted two BKT-type transitions for p = 6. 
A massive phase with a true long-range order appears 
below the lower transition at Tc2, while at intermediate 
temperatures Tc2 < T < Td a gapless phase with an 
algebraic quasi long-range order is realized. In our case 
the massive phase is represented by one of the states in 
Figs.[lja) and (b), while the gapless phase correspond to 
a state in Fig. [Ijc) with a power law decay of spin-spin 
correlations: 

iS!S^) ^ '-^^r^ ■ (13) 

The critical exponent rj continuously varies from ?;i = 1/4 
at T = Tci to = 1/9 at T = Tc2. The subsequent BKT 
transition related to the appearance of quasi-long-range 
order in the transverse components is expected to oc- 
cur at an independent transition temperature Tcs < Td . 
The expected sequence of finite-temperature phases is 
schematically shown in Fig. [3] with three BKT-type tran- 
sitions. A similar suggestion was made before for the tri- 
angular antiferromagnet with the exchange anisotropy;^ 
though no supporting numerical results were presented. 

To verify the outlined scenario in our case we have 
performed Monte Carlo simulations of the model ([T]) in 
zero magnetic field for the same value of the anisotropy 
parameter D/J = 1 as in Sec. II. Rhombic lattice clus- 
ters with periodic boundary conditions and with N — L"^ 
sites, i = 18 — 96, have been studied using the stan- 
dard Metropolis algorithm. Restricted motion of spins 
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FIG. 4; (Color online) Temperature dependence of the Binder 
cumulant for the antiferromagnetic order parameter mq for 
different cluster sizes. 



was implemented at low temperatures to keep the ac- 
ceptance rate around 50%. In order to improve further 
the performance of the MC algorithm, we have added a 
few microcanonical over-relaxation steps.— For mod- 
els without the single-ion term an over-relaxation move 
consists in a random rotation of a given spin about the 
local magnetic field. Such a step would not conserve the 
single- ion energy in ([T]). We choose, therefore, to reflect 
a spin with respect to the plane n-h, where n is the 
anisotropy axis and h is the local field. In total 2 • 10^ 
hybrid MC steps were used at each temperature and re- 
sults were further averaged over 20 different cooling runs, 
which both reduces measurement noise and provides an 
unbiased estimate of the statistical errors. 

The standard technique to locate a BKT transition is 
to measure the spin stiffness^ ^''^'''^^'■^^ which jumps from 
zero to the universal value Ps — 2TBKT/7r. However, in 
the case of an underlying discrete symmetry definition of 




FIG. 5; (Color online) Temperature dependence of the Binder 
cumulant for the uniform magnetization for different clus- 
ter sizes. 
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FIG. 6: (Color online) Scaling plot for the normalized sus- 
ceptibility versus the Binder cumulant in the vicinity of the 
upper transition Td. The indicated value of the exponent 
r) is used to achieve the best collapse of data from different 
clusters. 



FIG. 8: (Color online) Temperature dependence of the nor- 
malized susceptibility with 77 = 1/4 for different cluster sizes 
in the vicinity of the upper transition. The common crossing 
point yields T^i/J ^ 0.377. 



Ps becomes problematic. Therefore, we initially focus on 
the behavior of the Binder cumulant Ua — 
where A is the appropriate order parameter given by 
Eq. p4|) below. When correlations of the considered or- 
der parameter are critical the value of the Binder cumu- 
lant becomes size-independent. As a result, the curves 
Ul{T) measured for different cluster sizes L cross at the 
same point for a second-order transition, whereas for a 
BKT transition they merge once the correlation length 
is infinite. ■^'^ 

At every temperature we have separately measured 
even powers of different components of the order param- 
eter 

i,3 



for a — z and x, y and for q = Q, 0. Numerical results for 
z-components are presented in Figs. [3] and [51 which allow 
to locate approximately Td/J ~ 0.4 and Tc2/J ~ 0.2. 
For the second transition we use for illustration the uni- 
form magnetization instead of TOq. Nonzero values of 
unambiguously establish ^i-state in Fig. [IJa) as the 
low-temperature state with the broken Zg symmetry. In 
addition, this choice yields less noisy results. Still, sta- 
tistical errors are significant and the precise location of 
the transition point is difficult with this method. 

The renormalization group predictions^ for the expo- 
nent 77 in the vicinity of the two transitions can be, how- 
ever, tested without precise knowledge of the correspond- 
ing Tc-^^ In the critical regime the general scaling law^ 
reads as 

UUT) - fiL/0 and x - ^'""5(^0 , (15) 
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FIG. 7: (Color online) Scaling plot for the normalized sus- 
ceptibility versus the Binder cumulant in the vicinity of the 
second transition Tc2- The indicated value of the exponent 
rj is used to achieve the best collapse of data from different 
clusters. 
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FIG. 9: (Color online) Temperature dependence of the nor- 
malized susceptibility with 77 = 1/9 for different cluster sizes 
in the vicinity of the lower transition. The common crossing 
point yields Tc2/J ~ 0.198. 
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where x = ^^"^qa/^ is the generahzed susceptibihty, and 
^ is the correlation length. Hence, the plot of xl^'^~^ 
against Ul for the correct value of r\ should exhibit a 
collapse of numerical data for different cluster sizes onto 
a single curve. Figures [H] and [7] show the best fits around 
Tci and Tc2 respectively, which yield r\\ — 0.26 ±0.01 and 
r\2 = 0.12 ± 0.01. The obtained values are in a very good 
agreement with the prediction r\\ = 1/4 and 772 — 1/9.— 

Once the value of the critical exponent r\ is precisely es- 
tablished, one can use it to accurately estimate the tran- 
sition temperature from the finite-size scaling of suscep- 
tibility HID 4^1^ The curves x/^^"'' for different cluster 
sizes shown in Figs. |8] and [9] exhibit very tight crossing 
points giving us the following estimates for the transi- 
tion temperatures: T^/J = 0.377 ± 0.001 and T^ij J = 
0.198 ±0.001. 

The third BKT transition, which corresponds to a 
quasi-long-range ordering of transverse components, oc- 
curs at Tc3 < Tf-i- To precisely locate T^-i we measure the 
spin stiffness. The spin stiffness ps is defined as a general 
elasticity coefficient in response to a weak nonuniform 
twist of spins 8(\f^ (r) performed about a certain direction 
a in spin space. Generally, the spin stiffness is a fourth- 
rank tensor with the first pair of indexes running over 
the spin components and the second pair spanning over 
the gradient components in real space. In our case it is 
sufficient to consider only twists about the z-axis in the 
spin space, while all directions in the lattice plane are 
equivalent due to the six-fold rotational symmetry. This 
leaves us a single parameter: 



0.3 



6F^^ I d^^ 



(16) 



Choosing a twist with a uniform gradient along an ar- 
bitrary direction e in the lattice plane, one obtains in 
spherical coordinates 



Si • 



cos Ui cos ( 



sm t/,- sm ( 



: COs(<^i - (fj) 



with (fii = Lpi 



) e • r .; 



(17) 



Calculating the change of the free-energy up to the sec- 
ond order in a small 54> and normalizing result per unit 
area one obtains^Si^i^ 



2J2 



[Y.{S^Sy~SfS^)[e-{v,-v,)]]^) 



The first term in the above equation has been averaged 
over e = X and y directions. Numerical results from our 
MC simulations are presented in Fig. [TOl We determine 
crossing points of p^{T) with the straight line ps — 2T /n 
for each cluster size L and extrapolate them to i — > cxd 
according to T„oss{L) ^Tc3 + a/L. This yields the BKT 
transition at Tcs/ J — 0.168 ± 0.001 as illustrated in the 
inset of Fig. [TOl We have also determined the critical 
exponent 773 = 0.28 ± 0.03, which coincides within the 
error bars with the BKT value 77 = 1/4. 



T/J 




FIG. 10: (Color online) The spin stiffness of the easy-axis 
triangular antiferromagnet for different cluster size. Intersec- 
tion points with the line Ps/T — 2/tt are used to locate the 
transition temperature. The inset presents the interpolation 
of the crossing points to the thermodynamic limit. 



IV. SUMMARY 

We have studied a simple model of the Heisenberg 
triangular-lattice antiferromagnet with the single-ion 
anisotropy of the easy-axis type. Despite its simplic- 
ity such a 2D spin model exhibits a sequence of three 
BKT-type transitions illustrating nontrivial physical ef- 
fects which appear due to the competition between mag- 
netic anisotropy and geometrical frustration. The Monte 
Carlo simulations yield for D = J: T^i/J = 0.377, 
Tes/J = 0.198, and Tes/ J = 0.168. The two upper tran- 
sitions correspond to the breaking of the discrete Zg sym- 
metry, whereas the lowest one is the standard topological 
transition related to the proliferation of XY vortices. At 
Tc2 < T < Tci the longitudinal spin correlations have a 
power law decay with distance with a continuously vary- 
ing exponent 77. 

A remaining question is the fate of the intermediate 
critical phase at finite magnetic fields. An external field 
applied parallel to the anisotropy axis reduces the dis- 
crete symmetry from Zg to Z3. According to Jose et 
al-r^ the p = 3 clock model has no critical phase but ex- 
hibits instead a single transition into a normally ordered 
state, e.g., the phase in Fig. [IJa). It would be interesting 
to verify numerically the nature of this phase transition, 
which may be a critical one belonging to the three-state 
Potts model universality class^- or be of the first order 
due to a presence of the cubic term (jTO]) . 

The mean-field calculations find the partially ordered 
collinear phase. Fig. [Ijb), which appears to be unstable 
in 2D due to enhanced thermal fluctuations. Another in- 
teresting question for the future studies is whether the 
partially ordered state can be stabilized in layered tri- 
angular antiferromagnets. The thermal fluctuations are 
suppressed in this case by 3D effects, while the mean-field 
calculations predict stability of the partially disordered 
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phase up to Jj^ 0.6J. 
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